%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This is part of the set of files that accompany the article:       %
%Mankiw, N. Gregory and Ricardo Reis (2007) "Sticky Information in  %
%General Equilibrium," Journal of the European Economic Association,%
%forthcoming. See the appendix of the NBER or CEPR working paper    %
%versions for a detailed explanation of the algorithms.             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Please cite if you use the programs. I do not provide tech support.%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Last revised: August 30, 2006                                      %
%Written by: Ricardo Reis                                           %
%Input: the data in SIGEdata.mat                                    %
%Output: Bayesdraw#.mat with 50,000 draws form posterior            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%STEP 1: LOADING DATA AND MANIPULATING IT%%%%%
clear; clc; zxc=cputime;
load SIGEdata.mat; global T K; [T K]=size(data);
% Data matrix is Txk, where T is the number of observations and K is the
%number of variables in the model. First column is Dp, second is Dy, third
%is i, fourth isDwp, fifth is l.
global Y; Y=data(1,1:K)'; for i=2:T; Y=[Y;data(i,1:K)']; end;
% Y matrix is 5*Tx1 and has the observations stacked
global N; N=T;     %periods in MA representations, must be larger than T

%%%%STEP 2: PARAMETER INPUTS - CALIBRATED PARAMETERS %%%%%%%
beta=2/3;       %variable labor share
psi=4;          %wage elasticity of labor supply
theta=1;        %elasticity of intertemporal substitution
phiy=0.33;      %Taylor rule coefficient on output
phipi=1.24;     %Taylor rule coefficient on inflation
rho_m=0.9179;   %persistence of monetary policy shocks
sigma_m=0.0116; %stdev of shocks to money
rho_a=0.3496;   %persistence of tech growth shocks
sigma_a=0.0102; %stdev of shocks to tech
global pars_c;
pars_c=[beta psi theta NaN NaN phiy phipi rho_m sigma_m rho_a sigma_a NaN];
pars_c=[pars_c NaN NaN NaN NaN NaN NaN NaN NaN];

%%%%% STEP 3: BAYESIAN PRIORS %%%%%
%pars_e=[nu gamma rho_g sigma_g rho_nu sigma_nu rho_gam sigma_gam delta
%       omega lambda]
Pmean=[11 11 0.7 0.25^2 0.7 0.25^2 0.7 0.25^2 NaN NaN NaN]; %Prior means for par_e, note that for sigma it is the squared
Pvar= [10 10 0.05 0.25 0.05 0.25 0.05 0.25 NaN NaN NaN]; %Prior variances
%In Prior_par is a 2x11 vector with the 2 parameters characterising the
%distribution of each coefficient. For sigmas, it is v and s^2, for nu and
%gamma it is alpha and beta, for rho's it is alpha and beta
global Ppar
Ppar(2,1)=(Pmean(1)-1)/Pvar(1); Ppar(1,1)=Ppar(2,1)*(Pmean(1)-1);
Ppar(2,2)=(Pmean(2)-1)/Pvar(2); Ppar(1,2)=Ppar(2,2)*(Pmean(2)-1);
Ppar(1,3)=((1-Pmean(3))*Pmean(3)^2)/Pvar(3)-Pmean(3); Ppar(2,3)=Ppar(1,3)*(1/Pmean(3)-1);
Ppar(1,4)=4+(2*Pmean(4)^2)/(Pvar(4)); Ppar(2,4)=(Ppar(1,4)-2)*Pmean(4)/Ppar(1,4);
Ppar(1,5)=((1-Pmean(5))*Pmean(5)^2)/Pvar(5)-Pmean(5); Ppar(2,5)=Ppar(1,5)*(1/Pmean(5)-1);
Ppar(1,6)=4+(2*Pmean(6)^2)/(Pvar(6)); Ppar(2,6)=(Ppar(1,6)-2)*Pmean(6)/Ppar(1,6);
Ppar(1,7)=((1-Pmean(7))*Pmean(7)^2)/Pvar(7)-Pmean(7); Ppar(2,7)=Ppar(1,7)*(1/Pmean(7)-1);
Ppar(1,8)=4+(2*Pmean(8)^2)/(Pvar(8)); Ppar(2,8)=(Ppar(1,8)-2)*Pmean(8)/Ppar(1,8);
Ppar(1,9)=0; Ppar(1,10)=0; Ppar(1,11)=0; Ppar(2,9)=1; Ppar(2,10)=1; Ppar(2,11)=1;

%%%%% STEP 4: SIMULATION PARAMETERS %%%%%
%MLE estimates used for proposal function (with the sigma's)
x0=[34.068	4.196	0.93764	0.013929	0.62975	1.8194	0.66682	0.18657	0.18353	0.19516	0.70177]';
VARCOV = [0.99894	-0.073993	-0.00097616	-1.2784e-005	-0.0015209	0.035556	-3.9466e-005	-0.0050928	0.00033353	-0.00049702	0.0018762
-0.073993	0.39158	0.00095891	0.00039888	-0.00069195	-0.032885	-0.0083503	0.024084	-0.0068996	0.0013507	-0.0017569
-0.00097616	0.00095891	0.00043701	5.0562e-006	-0.00010796	0.001119	7.9993e-006	6.107e-005	-0.0002328	-0.00011831	9.0849e-005
-1.2784e-005	0.00039888	5.0562e-006	5.003e-006	-7.1753e-008	-8.7498e-005	-1.1179e-005	2.5167e-005	-1.1592e-006	-2.1911e-006	-7.7525e-006
-0.0015209	-0.00069195	-0.00010796	-7.1753e-008	0.00034887	0.0020227	0.00011452	3.3869e-005	0.00012301	-1.0153e-005	-0.00018504
0.035556	-0.032885	0.001119	-8.7498e-005	0.0020227	0.063531	-1.6795e-005	-0.0020255	0.0007714	-0.00085692	-0.00075835
-3.9466e-005	-0.0083503	7.9993e-006	-1.1179e-005	0.00011452	-1.6795e-005	0.0012026	5.5361e-006	-0.00011394	-8.83e-005	-4.8623e-005
-0.0050928	0.024084	6.107e-005	2.5167e-005	3.3869e-005	-0.0020255	5.5361e-006	0.0022174	-0.00054368	-3.0412e-005	-0.0002042
0.00033353	-0.0068996	-0.0002328	-1.1592e-006	0.00012301	0.0007714	-0.00011394	-0.00054368	0.00067406	-4.2356e-006	-5.8267e-005
-0.00049702	0.0013507	-0.00011831	-2.1911e-006	-1.0153e-005	-0.00085692	-8.83e-005	-3.0412e-005	-4.2356e-006	0.00012639	2.8861e-005
0.0018762	-0.0017569	9.0849e-005	-7.7525e-006	-0.00018504	-0.00075835	-4.8623e-005	-0.0002042	-5.8267e-005	2.8861e-005	0.00021824];
M=50000; burnin=30000; %Simulation parameters: M draws, burnin periods discarded
cc=0.75; kk=11;        %Drawing parameters: cc overdispersion, kk parameters
draw=[35.031	4.5894	0.92138	0.0128	0.62402	2.0682	0.6739	0.21896	0.17936	0.20829	0.72185];
        %Initial starting point. What changes across the 5 chains that I ran

 %%%%% STEP 5: DRAWING PARAMETERS FROM POSTERIOR %%%%%
logd=SIGEloglik(draw); prior=SIGEprior(draw); %useful starting values
draws=draw; accepts=1; logds=logd; %Storage of results
%Iterating over the Gibbs sampler steps
h = waitbar(0,'Going over Metropolis jumps iterations...');
for i=2:M
    [draw accept logd prior]=SIGEmetro(draw,(cc^2)*VARCOV,logd,prior);
    draws=[draws;draw]; logds=[logds;logd]; accepts=[accepts; accept];
    waitbar(i/M,h);
end
close(h); clear h draw accept logd temp prior; 
save Bayesdraw1.mat   %IMPORTANT: CHANGE THIS TO NOT OVERWRITE OLD RESULTS

%%%%% STEP 4: ANALYSE DRAWS %%%%%%
results = draws(burnin+1:M,:);
acceptsb = accepts(burnin+1:M,:);
avg = mean(results); vc = cov(results);
disp('Mean parameters and s.d.'); disp(avg); disp(std(results));
for i=1:kk
    subplot(kk,1,i), hist(results(:,i));
end
disp('rate_of_acceptance'); disp(mean(acceptsb))
disp('Routine took these seconds to run'); disp(cputime-zxc); clear zxc
